function generateMfuncMkernel4th(nameSteadyStateFile,allModelParams,y,x,xp,yp,strucOfArrays,position,nameOfFunction)
% This function generates an M function which evaluates all first, second, and third order 
% derivatives. The results are saved in the file nameOfFunction, for instance 
% nameOfFunction = 'numDerivM.m';


% We start deleting the old version of the file - if it exists
saveFile = [position,'\',nameOfFunction];
if exist(saveFile,'file') > 0
    delete(saveFile)
end

% Clears the command windue
%clc;
    
% The results are saved
%diary(saveFile);

text = ['function [M,Mxp,Myp,Mxpx,Mxpxp,Mxpy,Mxpyp,Mypx,Mypxp,Mypy,Mypyp] = ',nameOfFunction(1:end-2),'(params)'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Computing the steady ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = nameSteadyStateFile;
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'params = auxOut.params;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xss = auxOut.xssTrans;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'yss = auxOut.yssTrans;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Unfold params';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(allModelParams)
    text = [allModelParams{i}, '= params.',allModelParams{i},';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end 
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Current values in the model;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% x;';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(x(i)),' = xss(', num2str(i),');'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '% xp;';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(x(i)),'p',' = xss(', num2str(i),');'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% y';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
     text = [char(y(i)),' = yss(', num2str(i),');';];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '% yp;';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
     text = [char(y(i)),'p',' = yss(', num2str(i),');',];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Setting the dimension for the matrices at first order';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ny  = auxOut.ny;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx  = auxOut.nx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mxp   = zeros(1,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Myp   = zeros(1,ny);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Setting dimension for the matrices at second order';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mxpx   = zeros(nx,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mxpxp  = zeros(nx,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mxpy   = zeros(nx,ny);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mxpyp  = zeros(nx,ny);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');

text = 'Mypx   = zeros(ny,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mypxp  = zeros(ny,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mypy   = zeros(ny,ny);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Mypyp  = zeros(ny,ny);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');

names = fieldnames(strucOfArrays);
nn = size(names,1);
for i=1:nn
    DispSymMatrixMatlab_file(strucOfArrays.(names{i}),names{i},saveFile)
    text = ' ';
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
end

